转录组入门(6)： reads计数
往期教程
千万不要点阅读原文！
千万不要点阅读原文！
千万不要点阅读原文！
要求
实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。
需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来
对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。
看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等。
理论基础
在上篇的比对中，我们需要纠结是否真的需要比对，如果你只需要知道已知基因的表达情况，那么可以选择alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）。到了这一篇的基因（转录本）定量，需要考虑的因素就更加多了，以至于我不知道如何说清才能理清逻辑。
定量分为三个水平
基因水平(gene-level)
转录本水平(transcript-level)
外显子使用水平(exon-usage-level)。
在基因水平上，常用的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常用的HTSeq-count为例，这些工具要解决的问题就是根据read和基因位置的overlap判断这个read到底是谁家的孩子。值得注意的是不同工具对multimapping reads处理方式也是不同的，例如HTSeq-count就直接当它们不存在。而Qualimpa则是一人一份，平均分配。
对每个基因计数之后得到的count matrix再后续的分析中，要注意标准化的问题。如果你要比较同一个样本(within-sample)不同基因之间的表达情况，你就需要考虑到转录本长度，因为转录本越长，那么检测的片段也会更多，直接比较等于让小孩和大人进行赛跑。如果你是比较不同样本（across sample）同一个基因的表达情况，虽然不必在意转录本长度，但是你要考虑到测序深度（sequence depth)，毕竟测序深度越高，检测到的概率越大。除了这两个因素外，你还需要考虑GC%所导致的偏差，以及测序仪器的系统偏差。目前对read count标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之间的差异与换算方法已经有文章进行整理和吐槽了。但是，有一些下游分析的软件会要求是输入的count matrix是原始数据，未经标准化，比如说DESeq2，这个时候你需要注意你上一步所用软件会不会进行标准化。
在转录本水平上，一般常用工具为Cufflinks和它的继任者StringTie， eXpress。这些软件要处理的难题就时转录本亚型（isoforms）之间通常是有重叠的，当二代测序读长低于转录本长度时，如何进行区分？这些工具大多采用的都是expectation maximization（EM）。好在我们有三代测序。
上述软件都是alignment-based，目前许多alignment-free软件，如kallisto, silfish, salmon，能够省去比对这一步，直接得到read count，在运行效率上更高。不过最近一篇文献[1]指出这类方法在估计丰度时存在样本特异性和读长偏差。
在外显子使用水平上，其实和基因水平的统计类似。但是值得注意的是为了更好的计数，我们需要提供无重叠的外显子区域的gtf文件[2]。用于分析差异外显子使用的DEXSeq提供了一个Python脚本（dexseq_prepare_annotation.py）执行这个任务。
小结
计数分为三个水平： gene-level, transcript-level, exon-usage-level
标准化方法： FPKM RPKM TMM TPM
输出表达矩阵
在RNA-Seq分析中，每一个基因就是一个feature（特征？），而基因被认为是它的所有外显子的和集。在可变剪切分析中，可以单独把每个外显子当作一个feature。而在ChIP-Seq分析中，feature则是预先定义的结合域。但是确定一个read到底属于哪一个feature有时会非常棘手。因此HTSeq提供了三种模式，示意图见前一幅图
the union of all the sets S(i) for mode union. This mode is recommended for most use cases.
the intersection of all the sets S(i) for mode intersection-strict.
the intersection of all non-empty sets S(i) for mode intersection-nonempty.
基本用法非常的简单：
# 安装
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count用一个循环处理多个BAM文件(在/mnt/f/Data目录下)
mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log
done运行的时间会比较久，所以可以去了解不同参数的用法了，其中比较常用的为：
-f bam/sam： 指定输入文件格式，默认SAM
-r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序，默认是name
-s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的，所以需要判断到底来自于哪条链。如果选择了no， 那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上，第二个read则在另一条链上。
-a 最低质量， 剔除低于阈值的read
-m 模式 union（默认）, intersection-strict and intersection-nonempty。一般而言就用默认的，作者也是这样认为的。
-i id attribute: 在GTF文件的最后一栏里，会有这个基因的多个命名方式（如下）， RNA-Seq数据分析常用的是gene_id， 当然你可以写一个脚本替换成其他命名方式。
gene_id “ENSG00000223972.5_2”; transcript_id “ENST00000456328.2_1”; gene_type “transcribed_unprocessed_pseudogene”; gene_name “DDX11L1”; transcript_type “processed_transcript”; transcript_name “DDX11L1-002”; exon_number 2; exon_id “ENSE00003582793.1_1”; level 2;…
Jimmy的伏笔
我们这次分析是人类mRNA-Seq测序的结果，但是我们其实只下载了3个sra文件。一般而言RNA-Seq数据分析都至少要有2个重复，所以必须要有4个sra文件才行。我在仔细读完文章的方法这一段以后，发现他们有一批数据用的是其他课题组的： For 293 cells, the mRNA-seq results of the control samples include (1) those from the doxycycline-treated parental Flp-In T-REx 293 cells by us and (2) those from the doxycycline-treated control Flp-In T-REx 293 cells performed by another group unrelated to us (sample GSM1095127 in GSE44976)。 然后和Jimmy交流之后，他也承认自己只分析了小鼠的数据，而没有分析人类的数据。所以我们需要根据文章提供的线索下载另外一份数据，才能进行下一步的分析。
这个时候就有一个经常被问到的问题：不同来源的RNA-Seq数据能够直接比较吗？甚至说如果不同来源的RNA-seq数据的构建文库都不一样该如何比较?不同来源的RNA-Seq结果之间比较需要考虑 批次效应（batch effect) 的影响。
处理批次效应，根据我搜索的结果，是不能使用FPKM/RPKM，关于这个标准化的吐槽，我在biostars上找到了如下观点：
FPKM/RPKM 不是标准化的方法，它会引入文库特异的协变量
FPKM/RPKM has never been peer-reviewed, it has been introduced as an ad-hoc measure in a supplementary 没有同行评审
One of the authors of this paper states, that it should not be used because of faulty arithmetic 作者说算法有问题
All reviews so far have shown it to be an inferior scale for DE analysis of genes Length normalization is mostly dispensable imo in DE analysis because gene length is constant
有人建议使用一个Bioconductor包http://www.bioconductor.org/packages/devel/bioc/html/sva.html 我没有具体了解，有生之年去了解补充。
还有人引用了一篇文献 IVT-seq reveals extreme bias in RNA-sequencing 证明不同文库的RNA-Seq结果会存在很大差异。
结论： 可以问下原作者他们是如何处理数据的，居然有一个居然没有重复的分析也能过审。改用小鼠数据进行分析。或者使用无重复的分析方法，或者模拟一份数据出来，先把流程走完。
合并表达矩阵
HTSeq-count输出结果是一个个独立的文件，后续分析需要把多个文件合并成一个行为基因名，列为样本名，中间为count的行列式文件。肯定是不会用Excel手动处理（虽然可以写一个VBA脚本，但是数据量过大不好处理了），这里使用的Python写一个脚本。
基本逻辑：
读取文件
建立一个字典，如果key不在字典中，新增key和value，如果key在字典中，追加value。
输出
#!/usr/bin/python3
import sys
mydict = {}
for file in sys.argv[1:]:
for line in open(file, 'r'):
key,value = line.strip().split('\t')
if key in mydict:
mydict[key] = mydict[key] + '\t' + value
else:
mydict[key] = value
for key,value in mydict.items():
print(key + '\t' + value +'\n')这几行代码写了2个番茄钟，但是debug花了我一个番茄钟。问题出在str和list两种数据格式的混乱使用。还有一个bug： 由于词典是无序的，所以原本代表样本来源的第一行，会跑到其他行。
在论坛上找到一个更加简洁的代码（要求基因名顺序排列），用到paste,awk,printf这三个shell命令。paste *.txt | awk '{printf $1 "\t";for(i=2;i<=NF;i+=2) printf $i"\t";printf $i}'
保存为countCombiner.py，输入文件为count, 输出为标准输出，需要重定向。
简单分析
这一步需要用到R语言或者是Excel读取数据。
1.导入数据
options(stringsAsFactors = FALSE)# import data if sample are smallcontrol <- read.table("F:/Data/RNA-Seq/matrix/SRR3589956.count",
sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.count",
sep="\t",col.names = c("gene_id","rep2"))2.数据整合。gencode的注释文件中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜索到的，所以我就只保留ENSG00000105298这部分。
# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL3.总体情况, 大部分基因都为0，所以可以删掉节省体积。
summary(raw_count_filt)
control rep1 rep2
Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
Median : 0.0 Median : 0.0 Median : 0.0
Mean : 356.1 Mean : 370.3 Mean : 316.6
3rd Qu.: 15.0 3rd Qu.: 15.0 3rd Qu.: 10.0
Max. :161867.0 Max. :121902.0 Max. :105565.04.看几个具体基因
在EBI上搜索GAPDH找到ID为ENSG00000111640。
GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
gene_id control rep1 rep2
ENSG00000111640 ENSG00000111640.14_2 41857 53902 55302文章研究的AKAP95（ENSG00000105127）的表达量在KD中都降低了
> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
gene_id control rep1 rep2
ENSG00000105127 ENSG00000105127.8_2 1168 539 506下面的差异基因表达，让我想想，该如何收拾Jimmy挖的坑。
参考文献
[1] Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
[2] Detecting differential usage of exons from RNA-seq data
[3] RNA-seq Data Analysis-A Practical Approach(2015)